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PREDICTION  AND  ANALYSIS  OF  MATERIAL  RESPONSE  TO  IMPACT  AND  SHOCK  LOADING  USING 

A  SHARP-INTERFACE  EULERIAN  METHODOLOGY 
PI:  H.  S.  UDAYKUMAR 

DEPARTMENT  OF  MECHANICAL  AND  INDUSTRIAL  ENGINEERING 
UNIVERSITY  OF  IOWA,  IOWA  CITY,  IA-52242 


1.  Introduction 


We  have  developed  numerical  methods  [1-4]  and  a  computer  code  for  the  simulation  of 
multimaterial  interactions  that  result  from  high-speed  munition  impact.  This  is  an  important 
application  with  regard  to  the  Air  Force  with  implications  for  impact  and  penetration  of  targets, 
hazard  prevention  in  the  case  of  accidental  impact  on  explosives,  agent  defeat,  control  of  collateral 
damage,  control  of  atmospheric  spreading  of  hazardous  material  following  impact  etc.  The  direct 
beneficiary  agency  for  tins  work  is  the  AFRL-MNAC  at  Eglin  AFB,  FL.  This  work  was  initiated 
with  encouragement  from  Dr.  Kirk  Vanden,  chief  of  the  Computational  Mechanics  Branch  at 
AFRL-MNAC. 


Any  methodology  to  simulate  material  response  to  impact  and  shock  loading  should  be  able  to 
handle  the  following  physical  phenomena: 


1.  Large  deformations,  including  fragmentation  and  merger  of  the  materials. 

2.  Nonlinear  wave-propagation  and  the  development  of  shocks  in  materials  governed  by  rate- 
dependent  plasticity. 

3.  Accurate  elasto-plastic  deformation  of  the  interacting  objects  during  impact. 


A  fixed-grid  methodology  carries  the  advantage  of  allowing  for  arbitrary  deformation  while  avoiding 
problems  associated  with  moving  meshes.  This  implies  that  all  moving  bodies  will  be  embedded  in  a 
fixed  mesh  and  they  will  have  to  be  tracked  through  the  underlying  mesh.  Simultaneously,  the 
dynamics  of  such  moving  bodies,  i.e.  deformation,  collisions,  break-ups  etc.  will  have  to  be 
computed.  A  judicious  choice  of  methodology  to  solve  the  governing  equations  as  well  as  for 
tracking  the  moving  boundaries  has  to  be  devised  with  high  accuracy  and  robustness.  The 
characteristics  of  the  present  numerical  method  that  make  it  attractive  relative  to  other  methods 
employed  in  hydrocodes  for  high-speed  multimaterial  flows  are: 

1.  The  equations  governing  the  material  deformation  are  solved  in  an  Eulerian  setting  on  a 
fixed  Cartesian  mesh.  Well-developed  high-accuracy  shock  capturing  schemes  are  easily 
applied  to  compute  the  nonlinear  wave-propagation  phenomena  in  this  framework.  Here  the 
ENO  scheme  is  used  for  all  calculations.  Addition  of  problem-dependent  shock  viscosity  is 
not  called  for  since  adequate  dissipation  at  discontinuities  is  built  into  the  scheme. 

2.  The  mesh  remains  fixed  while  the  material  flows  through  the  mesh.  Thus,  issues  such  as 
mesh  deformation,  entangling,  catastrophic  mesh  distortion  in  regions  that  have  changed 
phase  from  solid  to  liquid  and  thus  have  lost  strength,  do  not  arise  within  the  current  fixed- 
grid  approach.  The  materials  can  fragment  or  collide  and/or  merge  without  affecting  the 
flow  calculations. 
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3.  The  interfaces  are  tracked  in  a  sharp  fashion.  They  are  not  smeared  over  the  mesh  as  in 
traditional  Eulerian  methods.  Thus  materials  can  approach  each  other  without  mixing  and  a 
mixture  formulation  is  not  required  in  treating  cells  with  multiple  interfaces  or  in  cells  that 
are  only  partially  filled.  The  exact  interface  location  is  known  at  all  times  due  to  the  level  set 
representation.  Boundary  conditions  and  jump  conditions  can  be  applied  at  the  sharp 
interfaces  at  both  free  surfaces  and  material-material  interfaces. 

4.  A  particle-level  set  method  is  used.  This  technique  is  shown  to  maintain  sharp  corners  of 
objects  without  excessive  smoothing  due  to  the  inherent  entropy  fix  in  the  advection 
scheme  used  to  evolve  the  narrow-band  grid-based  level  set.  Thus,  spurious  mass  loss 
effects  due  to  stretching  and  shearing  of  interfaces  that  plague  all  Eulerian  interface  tracking 
schemes  are  avoided  in  the  present  method.  No  difficulty  arises  in  treating  multiple 
boundaries.  These  are  simply  evolved  as  different  level  set  functions.  The  interfaces  can 
undergo  topological  changes  without  occasioning  difficulties  for  the  flow  solver. 

5.  Extension  to  3-dimensions  is  straightforward.  The  numerical  schemes  for  the  governing 
equations  are  obtained  by  field-by-field  decomposition  along  each  dimension  and  therefore 
addition  of  the  third  dimension  will  only  necessitate  discretization  of  the  equations  in  that 
direction.  Furthermore,  the  level  set  formulation  is  also  easily  extended  to  3D,  thus  allowing 
for  tracking  of  three-dimensional  objects  as  a  straightforward  extension  of  the  technique 
presented  in  this  work. 

2.  Objectives  and  Significance 

The  overall  project  goals  are  as  follows: 

1.  Advance  the  computational  techniques  for  simulation  of  high-speed  impact  on  fixed  Cartesian 
meshes  in  several  aspects,  viz.: 

(i)  Develop  a  local  refinement  strategy  to  improve  the  accuracy  and  efficiency  of  the 
method.  Cartesian  methods  suffer  from  lack  of  flexibility  in  placing  grid  points  in 
locally  refined  regions.  Thus,  when  the  mesh  needs  to  be  refined  to  improve 
accuracy  or  capture  small-scale  features,  the  grid  density  increases  globally.  To 
overcome  this  drawback,  local  refinement  techniques  are  being  developed. 

(ii)  Implement  damage  and  materials  models  to  be  able  to  simulate  material  perforation 
and  fragmentation.  Rate-dependent  constitutive  modeling  such  as  in  the  Johnson- 
Cook  model  to  model  visco-plastic  deformation  has  been  implemented. 

(iii)  Extend  the  method  to  3-dimensions.  This  is  relatively  straightforward  to  do  within 
the  framework  of  level-set-based  interface  representation.  To  render  3D  simulations 
feasible  parallelization  of  the  code  is  imperative.  The  code  is  currently  operational  in 
parallel  mode  and  further  refinements  are  being  performed  to  optimize  performance. 

2.  Apply  the  method  to  problems  of  interest  to  collaborators  at  AFRL-MNAC  (Eglin  AFB,  FL).  The 
specific  types  of  problems  to  be  investigated  include: 

(i)  The  perforation  and  damage  of  metals  under  high  strain-rate  conditions,  such  as 
under  the  effects  of  shock  loading,  impact  and  detonation  wave  loading.  Of  particular 
interest  are  the  micro-scale  damage  mechanisms,  including  void  formation,  growth 
and  coalescence  leading  to  spall  in  ductile  materials. 

(ii)  Collapse  of  voids  in  energetic  materials  due  to  imposition  of  impact/ shock  loading. 
This  problem  is  of  interest  in  high-energy  density  explosives  where  accidental  impact 
can  initiate  detonation  and  catastrophic  explosion. 


3.  Work  accomplished  in  this  project 

3.1  Simulations  of  high-speed  multimaterial  Interactions 


For  2-dimensional  problems  the  hybrid  particle  level  set  method  has  been  used  to  track  boundaries 
with  sharp  comers  that  are  carried  without  deterioration  through  the  large  deformations  of  the 
materials.  The  details  of  the  method  have  been  presented  in  two  journal  papers  [1-2].  Benchmark 
calculations  for  the  multi-dimensional  case  including  axisymmetric  Taylor  bar  impact  and 
penetration  of  a  Tungsten  rod  into  steel  plate  have  shown  excellent  agreement  with  moving  finite 
element  solutions.  Qualitative  agreement  with  theory  is  shown  for  void  collapsing  process  in  an 
impacted  material  containing  a  spherical  void.  The  method  has  thus  been  shown  to  be  suitable  for 
applications  involving  high-velocity,  multimaterial  impacts  leading  to  large  strain-rates,  nonlinear 
elasto-plastic  waves  and  topological  changes. 

An  example  of  the  capability  [2]  of  the  method  is  shown  in  Figure  1 .  The  validation  of  our  method 
for  two  deformable  objects  with  different  material  properties  (a  case  requiring  2  different  level  sets) 
is  carried  out  using  a  slender  tungsten  heavy  alloy  (WHA)  rod  projectile  penetrating  an  initially 
planar  target  made  of  a  steel  plate  with  a  velocity  of  1250  m/s.  A  Johnson-Cook  material  model  is 
used  and  the  corresponding  strength  parameters  for  both  materials.  Note  that  friction  between  the 
two  impacting  surfaces  is  neglected  in  these  calculations.  The  evolution  of  equivalent  plastic  strain  is 
shown  in  Figs.  l(a-c)  for  the  extended  domain  with  a  160x688  mesh.  The  maximum  equivalent 
plastic  strain  is  found  to  be  around  4.5,  occurring  mostly  near  the  impact  surfaces.  The  values  of 
equivalent  plastic  strain  are  higher  in  the  WHA  material  compared  to  those  in  the  steel  material. 
The  plastic  strains  obtained  by  Camacho  and  Ortiz  [7]  using  Lagrangian  finite  element  method  with 
an  adaptive  mesh  agree  very  well  with  the  present  results,  both  in  terms  of  the  magnitude  and 
distributions  of  the  plastic  strains.  In  particular,  a  trough  in  the  plastic  strain  distribution  is  noticed 
in  both  our  results  as  well  as  those  of  Camacho  and  Ortiz  and  occurs  near  the  bottom  surface  in  the 
steel  plate  at  the  symmetry  axis,  as  seen  in  Fig.  1(b).  The  ejection  length  of  the  WHA  material  is 
higher  in  the  Camacho  and  Ortiz  calculations  when  compared  to  our  results.  However,  the 
resolution  of  the  ejected  region  afforded  by  the  mesh  used  in  the  present  calculations  is  too  low, 
with  just  3  mesh  points  across  the  vertically  oriented  trails  of  the  ejecta.  The  grid  refinement  study 
performed  above  indicates  that  as  the  mesh  is  refined  further  the  length  of  the  ejecta  will  increase. 
As  shown  below,  at  the  current  mesh  resolution,  the  overall  penetration  characteristics  and  material 
deformation  are  adequately  predicted.  Fig.  2  shows  the  projectile  nose  and  tail  trajectories  as  a 
function  of  time,  for  the  extended  domain  case,  and  is  compared  with  the  superposed  results  from 
experimental  data  and  from  Camacho  and  Ortiz  [7].  Also  plotted  are  the  original  rear  and  impact 
surfaces.  Our  results  show  reasonable  agreement  with  those  of  experiment  and  Camacho  and  Ortiz. 
The  tail  trajectory  is  in  much  better  agreement  as  its  surface  experiences  less  extreme  conditions 
during  impact  and  penetration.  The  present  calculation  predicts  the  penetration  depth  in  good 
agreement  with  experiments.  Despite  the  marginal  resolution  of  the  ejected  trails,  the  overall 
penetration  and  deformation  behavior  is  predicted  in  good  accord  with  the  adaptive  finite  element 
simulations  of  Camacho  and  Ortiz. 


3.2  Investigation  of  thermomechanics  of  void  collapse  in  an  energetic  material 


In  the  recently  completed  PhD  work  of  Linhbao  Tran,  the  methodology  was  advanced  to  include 
the  thermomechanics  of  void  collapse  in  an  energetic  material  that  is  subject  to  shock  loading.  The 
Eulerian,  sharp  interface,  fixed  Cartesian  grid  method  was  applied  to  study  hot  spot  formation  in  an 
energetic  material  (HMX)  under  various  loading  conditions  [3-4].  The  mass,  momentum,  and  energy 
equations  were  solved  along  with  evolution  equations  for  deviatoric  stresses  and  equivalent  plastic 
strain.  These  equations  were  cast  in  Eulerian  conservation  law  form.  Pressure  was  obtained  from 
the  Mie-Griineisen  equation  of  state.  The  material  was  modeled  as  a  viscoplastic  solid.  High-order 
accurate  ENO  shock-capturing  schemes  along  with  a  particle  level  set  technique  were  used  to  evolve 
sharp  immersed  boundaries.  The  details  of  void  collapse  under  shock-loading  and  the  coupling  of 
the  void  dynamics  to  the  energy  release  in  the  solid  material  were  analyzed  and  are  presented  in  the 
papers  [3-4].  The  thermo-mechanical  response  of  solid  phase  was  combined  with  the  vapor-phase 
compression  and  chemical  reactions  to  study  their  effects  on  the  energy  deposition  mechanisms. 
The  formulation  and  numerical  treatment  of  the  coupled  solid-gas  interactions,  including  heat 
release  due  to  chemical  reactions  was  also  successfully  tackled.  The  effects  of  loading  intensity  and 
void  size  were  studied.  The  results  show  that  for  the  micron-size  voids  under  consideration 
significant  gas  phase  chemical  reactions  occur.  However,  their  influence  on  the  void  collapse  itself  is 
minimal.  Figure  3  shows  a  sample  result  for  the  collapse  of  a  20  micron  diameter  void  when  the 
lower  boundary  of  the  domain  is  subject  to  shock  loading  with  imposed  particle  velocity  of 
Up=100m/s.  When  the  shock  collapses  there  is  a  hot  spot  created  where  significant  chemical 
reactions  occur  leading  to  the  breakdown  of  HMX  into  monomer  and  other  reactive  components. 
The  point  of  collapse  is  also  a  region  of  high  temperature  and  pressure  and  a  compression  wave 
emanates  from  the  collapse  point  to  the  surrounding  material.  Results  [3-4]  indicate  that  for  micron- 
size  voids  the  collapse  time  is  too  short  to  set  up  thermal  runaway  at  the  hot  spot  prior  to  void 
collapse. 

3.3  Extension  to  multimaterial  interactions  in  low-speed  flows 

As  an  extension  of  the  method  development  research  that  is  being  undertaken  in  this  project, 
advances  have  resulted  in  the  development  of  the  level-set-based  technique  for  handling  of  low- 
speed  (incompressible  flow  regime)  fluid-solid  interaction  problems.  The  attractive  feature  of  the 
current  research  is  the  generality  that  it  affords  in  handling  a  wide  range  of  multimaterial  problems 
with  moving  boundaries  in  3-dimensions.  In  recently  submitted  papers  [5-6]  we  have  demonstrated 
the  capability  of  the  method  to  tackle  fluid-structure  interaction  problems  [5]  and  droplet-surface 
interaction  problems  [6].  In  these  papers  the  ease  of  implementation  of  the  technique  in  3- 
dimensions  for  such  problems  has  also  been  demonstrated.  These  applications  were  made  possible 
by  the  general  idea  of  using  level-sets  to  track  boundaries  and  the  ability  to  handle  collisions  of 
interfaces.  The  central  theme  of  the  current  research  effort,  i.e.  to  develop  a  sharp-interface 
discretization  strategy  has  been  used  in  these  papers.  Validation  of  the  numerical  results  has  been 
performed  in  each  case  with  experimental  data  and  other  numerical/analytical  solutions. 

A  fluid-structure  interaction  problem,  that  of  a  sphere  oscillating  in  a  stratified  fluid,  is  shown  as  an 
example  in  Figure  4.  For  a  specific  set  of  sphere  oscillation  frequency  and  Reynolds  number  of  the 
oncoming  flow,  the  vortical  structures  visualized  using  the  A-method  are  plotted  in  Figure  4.  The 
cylindrical  vortices  apparent  in  the  x-y  view  shown  in  Figure  4(b)  illustrate  that  the  vortex  shedding 
is  self  similar  in  z-direction.  This  demonstrates  the  quasi-two-dimensionality  imposed  due  to 
stratification  and  is  in  perfect  agreement  with  the  experiments  on  this  problem  [8].  The  x-z  view  of 
vortex  shedding  shown  in  Figure  4(c)  supports  the  above  observation.  The  Strouhal  number 


calculated  from  the  probe  velocity  variation  with  time  is  0.35  which  indicates  that  the  shedding 
under  these  conditions  locks  on  with  the  oscillation  frequency  of  the  sphere.  Figure  4  also  illustrates 
that  the  vortices  are  being  shed  alternately  from  either  side  of  the  sphere.  Hence  this  flow  regime  is 
named  “lock-on  alternate  single  vortex”  regime  [8].  The  flow  structures  viewed  from  the  various 
perspectives  are  visually  identical  to  those  obtained  in  the  experiments  [8].  In  summary,  flow  around 
an  oscillating  sphere  in  the  presence  of  density  stratification  has  been  simulated  for  a  selected  set  of 
parameters.  The  flow  regimes,  vortical  patterns  and  shedding  frequencies  predicted  by  the  current 
technique  for  these  test  points  are  in  agreement  with  the  experimental  results  reported  by  Lin  & 
coworkers  and  thereby  provide  validation  for  the  present  methodology  for  three-dimensional 
incompressible  flows  involving  moving  immersed  objects. 

A  sample  result  of  multimaterial  incompressible  flows  is  shown  in  Figure  5  (see  [6]  for  extensive 
validation  of  the  method  for  such  applications).  The  calculation  carried  out  is  for  a  water  droplet 
suspended  in  air  which  falls  on  a  stationary  solid  cylindrical  surface.  The  following  dimensionless 
parameters  apply:  Reynolds  number  Rtf=10,  Weber  number  We—  1 0,  and  Froude  number  Fr=  1.  The 
liquid  droplet  (density  of  the  liquid  is  100  times  that  of  the  surrounding  air)  starts  to  move  down  due 
to  gravity,  impacts  on  the  curved  surface,  spreads  and  drips  down,  as  observed  in  Figure  5. 
Eventually,  this  droplet  undergoes  severe  deformation  following  which  the  droplet  breaks  and  forms 
two  smaller  drops  that  fall  under  gravity.  The  falling  film  of  fluid  then  retracts  due  to  capillary  forces 
and  accumulates  to  form  a  second  drop,  which  in  turn  is  pulled  down  due  to  gravity  and  breaks 
away  from  the  main  drop.  This  problem  shows  the  capability  of  the  method  to  track  interactions 
and  severe  deformations  of  boundaries  in  low-speed  flows. 

4.  Work  in  Progress  and  extensions 

The  capability  to  compute  the  flows  in  3D  has  been  developed  in  the  framework  of  a  fixed  Cartesian 
mesh  and  level-set  representation  of  the  interfaces.  In  order  to  make  3D  computations  feasible,  the 
code  has  been  parallelized  using  MPI  (Message  Passing  Interface).  Local  mesh  refinement  using  an 
octree  structure  to  subdivide  the  mesh  points  has  also  incorporated  into  the  code.  The  use  of  local 
mesh  refinement  produces  issues  that  conflict  with  parallel  efficiency  and  thus  care  needs  to  be 
exercised  in  coordinating  the  parallel  implementation  with  mesh  refinement  strategy.  Furthermore 
the  behavior  of  waves  at  fine-coarse  mesh  interfaces  needs  to  be  carefully  treated  or  spurious  wave 
reflections  or  damping  will  result.  These  challenging  issues  are  being  addressed  in  current  work. 
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Figure  1.  Simulation  of 
penetration  of  a  Tungsten  rod 
into  a  steel  plate.  Impact 
velocity  is  1.25  km/s.  The 
shapes  of  the  rod  and  plate  are 
shown  at  three  instants  of 
time.  Contours  of  plastic 
strain  are  also  shown. 
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Figure  2.  Comparison  of  the 
positions  of  the  tail  and  head 
of  the  Tungsten  rod  shown  in 
Figure  1  as  obtained  from  the 
current  calculation  and  the 
numerical  (FEM)  benchmark 
[8]  and  experimental  data. 
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Figure  3.  (a)  Evolution  of  void  collapse  process  for  void  diameter  of  20  microns,  impact 
velocity  of  100  m/s  imposed  on  the  lower  surface  of  the  computational  domain,  shock  rise 
time  of  10  ns.  (b)  Temperature  field  at  point  of  collapse,  (c)  Decomposed  HMX  species 
concentration  field  at  point  of  void  collapse. 
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Figure  4.  Results  for  3- 
dimensional  fluid-solid 

interaction  problem.  A  sphere 
oscillates  in  the  streamwise 
direction  with  an  oncoming  flow 
in  a  stratified  fluid  (stratification 
is  along  z-direction).  The  figures 
show  vertical  structures  at  a 
Reynolds  number  of  190  for  a 
specific  oscillation  frequency  of 
the  sphere.  See  [3]  for  details, 
(a)  Oblique  view,  (b)  x-y  view, 
(c)  x-z  view. 
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Figure  5.  Shapes  of  a  droplet  impacting  on  a  cylinder  (Re=10,  We=333).  The 
formation  of  a  filament  on  each  side  is  seen  along  with  the  detachment  of  a 
pendant  drop  and  formation  of  a  secondary  pendant  drop  after  retraction  of  the 
remaining  filament. 


